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Abstract 

' This paper proposes a new methodological framework within which the heat conductance in ID 



o 
o 



X 



lattices can be studied. The total process of heat conductance is separated into two parts where 
the first one is the equilibrium process at equal temperatures T of both ends and the second one 
— non-equilibrium with the temperature AT of one end and zero temperature of the other. This 
approach allows significant decrease of computational time at AT — t- 0. The threshold temperature 



Q 

U, 

. Tthr is found which scales Tthr(-^) ~ ^ with the lattice size N and by convention separates two 



mechanisms of heat conductance: phonon mechanism dominates at T < Tthr and the soliton 
contribution increases with temperature at T > Tthr- Solitons and breathers are directly visualized 
in numerical experiments. The problem of heat conductance in non-linear lattices in the limit 



> 

, AT — )• can be reduced to the heat conductance of harmonic lattice with time-dependent stochastic 



rigidities determined by the equilibrium process at temperature T. The detailed analysis is done 
for the /3-FPU lattice though main results are valid for one-dimensional lattices with arbitrary 
potentials. 
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I. INTRODUCTION 



The problem of heat conductance in low dimensional systems attracts much attention 
in last decades (see review and is motivated by the discovery of quasi-one-dimensional 
(nanotubes, nanowires, etc.) and two-dimensional (graphen, graphan, etc.) systems. 

The modern theory of heat conductance was initiated by the celebrated preprint of 
E. Fermi, J. Pasta and S. Ulam though the primary aim was "0/ establishing, exper- 
imentally, the rate of approaching to the equipartition of energy among the various degrees 
of freedom" . Subsequent investigations demonstrated wide area of consequences in many 
physical and mathematical phenomena (see reviews in special issues of journals CHAOS |3| 
and Lecture Notes in Physics j4| devoted to the 50th anniversary of the FPU preprint). 

The dynamical properties of nonlinear systems in microcanonical ensemble (total energy 



E = const) were thoroughly ana 
and to get exact results (soliton 



yzed in most papers. It allows to investigate the dynamics 
5-7| and breather js-jj] solutions), to analyze regular and 



stochastic regimes and to find the corresponding thresholds. The 



'PU preprint also initiated 
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the investigations in the field of "experimental mathematics" 

About ten decades ago P. Debye argued that the nonlinearity can be responsible for the 
finite value of heat conductance in insulating materials 1^. But modern analysis shows 
that it is not always the case. There are many examples where the coefficient of heat 
conductance x diverges with the increasing of the system size L as x oc where a > 0, 
and X —i- 00 in the thermodynamic limit {L 00). Most of momentum conserving one- 
dimensional nonlinear lattices with various types of nearest-neighbor interactions have this 
unusual property (see, e.g.^, [is], Q ). Moreover, some other systems, - two- 17;20| and 



three-dimensional lattices 



21 



22j, polyethylene chains 23|, carbon nanotubes 24N28| have 



analogous property - diverging heat conductance with the increasing size of the system. 

There were some conjectures explaining the anomalous heat conductance. Generally 
speaking, whenever the equilibrium dynamics of a lattice can be decomposed into that of 
independent "modes" or quasi-particles, the system is expected to behave as an ideal thermal 



conductor 



29| . Thereby, the existence of stable nonlinear excitations is expected to yield 



ballistic rather than diffusive transport. At low temperatures normal modes are phonons. 
At higher temperatures noninteracting "gas" of solitons starts to play more significant ro 
and M. Toda was the first, who suggested the possibility of heat transport by solitons 30 1. 
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Though analytical expressions for solitons can be derived only for few continuum mode 



described by partial differential equations, Friesecke and Pego in a series of recent papers 
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341 ] made a detailed study of the existence and stability of solitary wave solutions on discrete 



lattices with the Hamiltonian H = + u{yi) , where yi = Xi — Pi = Xj. It has been 

proven that the systems with this Hamiltonian and with the following generic properties of 
nearest-neighbor interactions: u'{0) = 0; ^"(O) > 0; u"'{0) ^ has a family of solitary wave 
solutions which in the small amplitude, long-wavelength limit have a profile close to that of 
the KdV soliton. It was also shown [35] that these solutions are asymptotically stable. Thus 
most acceptable point of view on the origin of anomalous heat conductance in nonlinear 
lattices is as follows: phonons are responsible for heat conductance at low temperatures, 
and at high temperatures - solitons [36|, l37 |. 

A set of generic properties were found in a series of papers in investigation of dynamics of 
nonlinear lattices, starting from the celebrated preprint of FPU And one is an existence 
of stochasticity thresholds. The weak stochasticity threshold is characterized by a specific 
energy S below the which the trajectory in the phase space is almost regular (with near 
zero Lyapunov exponents) and only small part of normal modes is excited (it is just the case 
observed and analyzed by FPU). The strong stochasticity threshold corresponds to the value 
of £ above which energ y equip artition between normal modes is established, and Lyapunov 



exponents are positive 
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A major part of results was obtained using microcanonical ensemble for isolated systems. 
Physically more justified is the usage of canonical ensemble where temperature is kept (on 
average) constant by some or other type of heat baths. If the constant temperature is 
maintained by the Langevin sources (random forces with viscous friction), then from the 
Fokker-Planck equation the equilibrium Gibbs distribution immediately follows. 

If one starts calculations from arbitrary initial conditions in canonical ensemble then some 
time is necessary to achieve the state of thermodynamic equilibrium. And this stage is not a 
trivial one 46|] . Firstly and surprisingly, kinetic and potential energy can relax to equilibrium 
with different rates; secondly, obeying the Maxwell velocity distribution function is not the 
sufficient condition of achievement the equilibrium. And the critical stage of achievement 
the equilibrium (energy equipartition between normal modes) is the excitation of the most 
longwave normal mode. Characteristic times r of achievement the equilibrium can cover 
very wide range. For instance, there are well localized excitations in the harmonic lattice 
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with random masses or, equiyalently, random interparticle potentials (Anderson localization 
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48[) where r > 10 



300 



491]. And this phenomenon is explained by very weak interaction 



of localized excitations, if they are centered near the lattice center, with the heat reservoir 
located at the lattice ends. 

If to return back to the problem of heat conductance, then one meets rather confusing 
experimental and numerical results, e.g. exponent a in the dependence x oc depends on 
the model under consideration, types of boundary conditions, used thermostat (Langevin 



or Nose-Hoover 



50 



51|). and also on temperature. For instance, temperature dependence 



of heat conductance in carbon nanotubes decreases as x ~ 1/T at T > 10 K [52|]; exper- 
imentally is found [25] that x also decreases with the growth of temperature. Different 
temperature dependencies x vs. T were found in ID nonlinear lattices. For /3-FPU lattice: 

n 

X ~ N°'T ^ at T < 0.1 and x ~ JSl^T^I'^ at T > 50 [53] what is usually observed in insulating 
crystals. For the interparticle harmonic potentials and on-site potentials (e.g. Klein-Gordon 
chains) x ~ T^^'^^, i.e. heat conductance decreases with the growth of temperature 54 1. 
One more problem is the calculation of heat conductance at small temperature gradients. 
Usually these calculations are very time consuming because of great fluctuations of heat 
current and statistical averaging over large number of MD trajectories is necessary. 

The paper organized as follows: in Section II we introduce new method for the calculation 
of the heat conductance which significantly decreases the computation time and diminishes 
the standard error. The method is based on the separation of the total process of heat 
conductance into two contributions: equilibrium and non-equilibrium and the latter one 
is responsible for the energy transfer. In the next Section we found that some quadratic 
mean values do not exhibit the expected tendency to reach zero values as the temperature 
difference AT — ?■ 0. The threshold temperature Tthr, separating two regimes, - damped 
and undamped, is revealed. And the dependencies of Tthr on temperature T and lattice 
length are found. Some modifications of the calculation of heat conductance in the limit 
AT — > are introduced in Section V. Direct evidences of the solitons contribution to the 
heat conductance are given in the next Section. /3-FPU lattice is considered as an example. 
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II. HEAT CONDUCTANCE IN THE /3-FPU LATTICE 



We consider the one-dimensional lattice of oscillators with the interaction of nearest 
neighbors 

U = '^u{yi), yi = Xi-Xi_i (1) 

i 

and the /3-FPU potential u{y) = + jy^ (usually we put m = /3 = 1). 

Nonequilibrium conditions are necessary for the heat transport simulation. The most 
abundant method is the placement of the lattice into the heat bath with different tem- 
peratures of left T+ and right T_ ends (T+ > T_). Different types of heat reservoirs are 
thoroughly analyzed in We utilize the Langevin forces acting on the left = — 7x1 
and right F^ = —'-/xn oscillators. {^±} are independent Wiener processes with zero mean 
and (^±(^1) '^±(^2)) = 27T±(5(ti — ^2)- AT = (T+ — T_) is the temperature difference. The 
generalized Langevin dynamics with a memory kernel and colored noises is also suggested 



55( 1 to correctly account for the effect of the heat baths. 

The following set of stochastic differential equations (SDEs) 

dU 

^i = —Q^.+ + ^i^F- (2) 

are to be solved to find the heat flux J. Then from the Fourier low J = — x VT the coefficient 
of heat conductance is 

X = NJ/AT, (3) 
and the problem is to find the heat current J. The local heat flux (from ith to {i + l)th 



oscillator) is defined [56| by 



Ji-^i+i — {Fi^i+i Xi+i) ; = —U'{xi^i — Xi), (4) 

where -Fj^j+i is a shorthand notation for the force exerted by the ith on the {i + l)th 
oscillator and (...) is the time averaged. The total heat flux J can be found as the mean 
value J=(iV-l)-i^f-ij,^,+i. 

A. Equilibrium and non-equilibrium contributions to the heat conductance 

If r_ 7^ then the process of heat conductance can be formally separated into two 
parts: the first one - equilibrium process with equal temperatures T_ of both lattice ends; 
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FIG. 1: Schematic representation of the sphtting of the total process x(t) into equihbrium ^^{t) 
and non-equiUbrium x^(t) ones. 

and second - nonequilibrium process with temperature AT of the left lattice end and zero 
temperature of the right end (see Fig. [T]) (by 'process' we hereafter assume for brevity the 
solution x(t) = Xi(t),X2{t), . . . ,XN{t); v(t) = Vi(t),V2{t), . . . ,VN{t) of the corresponding 
SDKs). 

Namely the second process defines the heat transport realized against the background 
of the equilibrium process. Once we utilize this approach then the Langevin forces in ([2]) 
can be written as {^+} = {^°} + {^^} and {^-} = {C,'^} for the left and right lattice ends, 
correspondingly; superscripts '0' and '1' refer to equilibrium and nonequilibrium processes. 
Then the total dynamical process x(t) can be represented as the sum of two processes 



x(t) =x°(t) + xi(t), 



(5) 



where x''(t) is the equilibrium (Gibbs's) process at temperature T_, and x^(t) - nonequilib- 
rium, responsible for the energy transport, process. Then the Langevin dynamics is 



dx 

dU dU° 

dx, 



+ 5a(e°-7i?) + 5^7v(e°-7i°^), 



(6) 



dxi 



+ 6iii^^ -^x\) + 6iNi-lxlj), (7) 

and the sum of equations and ([7]) is virtually identical to the parent equation 
dg). Random values {^°} and {^^} obey the identities (^°(ti)^°(ti)) = 2jT_6{ti - tg) 
and {^^(ti)^^{ti)) = 2jAT6(ti — ^2); U'^ is the total energy ([T]) where the arguments 
xi{t),X2(t), . . . ,XAr(t) of the total process are substituted to the coordinates of the equilib- 
rium process x5(t), ^^(t), . . . , a;5^(t). Expression in the square brackets in ([7]) is the difference 



of forces acting on the ith particle from the total process x(t) and equilibrium process x°(t). 
It is significant to note that this force is the random value, and the process x^(t) (heat 
transport) is realized in the lattice with time- dependent random potentials. The problem of 
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heat conductance in the random time-independent potentials was analyzed in 

Equation describes the system embedded in the heat reservoir at temperature T_. 
And x°(t) is the stationary equilibrium process described by the canonical Gibbs distribution 
(equilibrium thermodynamics of the /9-FPU lattice in the canonical ensemble was considered 



m 



Process x^(t) is responsible for the heat transport and the Wiener's process {'C^(t)} on 
the left lattice end defines small temperature AT. Right lattice end has zero temperature. 
An expression for the local heat flux is 

Ji^i+i = (-Fi_^i+i(x) Xi+i - Fi_,i+i(x°) , (8) 

and the equilibrium process x° does not transfer energy: ^Fj^j+i(x°) = 0. 

One of the goals of the present paper is the calculation of heat conductance at small 
temperature gradients. Usually these calculations are realized by solving SDEs ([2]) and are 
very time consuming because of great fluctuations of heat current (below we show that the 
time of computation increases cx (AT)"' if the accuracy of calculations is predetermined). 

The comparison of two approaches (solving of standard SDEs and (I6])-(I7!)) is shown 
in Fig. [2] and results coincide with very good accuracy. Note, that most of results in this 
paper are presented for the number of oscillators = 5 in the lattice. It may appear that 
this value is too small. For instance, the best estimate so far required simulations of up 
> 10'^ particles and > 10^ integration steps plus ensemble averaging But our results 
are aimed at founding some basic issues where number of particles is less essential. Lattices 
with larger number of oscillators were tested where necessary. 

The dependence of heat conductance on the particles number A^ is shown in Fig. [3] at 
two value of temperature T_. Inharmonicity becomes negligible in the limit T_ — )■ and the 
analytical solution of the heat conductance for the harmonic lattice is given in js^. 

There should be solved twice as large SDEs (I6])-(I7]) in suggested approach as that in 
standard scheme ([2]), and this the price which is paid for the facility with using small 
temperature gradients. As one would expect, the accuracy of the suggested approach is 
higher (provided that all computational terms and conditions are identical). The comparison 



1.5 




3 T 4 



FIG. 2: Coefficient of heat conductance vs. temperature for the lattice of = 5 oscillators. 
Filled circles: solution of standard SDEs ([2|); empty circles: SDEs ©-([Tl). Averaging over 100 MD 
trajectories 10^ time units (t.u.) each. T_ = 0.2, AT = O.Oir_, 7 = 1. Triangle up at T = is 
the exact value in the harmonic approximation (/3 = 0). 

of accuracies is given in Appendix |X] 

III. STRANGE BEHAVIOR OF PROCESS x.^{t) AT HIGH TEMPERATURES 

Usually the temperature difference AT ~ (0.01 — 0.1)T_ at the lattice ends is an ap- 
propriate choice. Then the Fourier law J oc AT (at fixed A^) is valid with good accuracy. 
Actually, the corrections to the heat current are of the order (AT)^ as the current is the 
odd function of the temperature difference, and this ensures the reasonable accuracy of the 
linear approximation. 

Now we concentrate our efforts on the elucidating the heat conductance dependence 
via temperature of the background process x°(t). Langevin forces {^^}, which provide 
temperature AT, are of the order ~ V AT (as (^^(ti)^^(t2)) ~ AT). And one can expect 
that process x^(t) should have the same order x^(t) ~ V AT because equation ([7]) becomes 
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FIG. 3: Coefficient of heat conductance for the /3-FPU lattice for N = 7 — 150 oscillators. Squares: 
T_ = 1, circles: T_ = 0.1. Filled symbols - results obtained by the solution of standard SDEs ([2]), 
empty symbols - SDEs dS])-©. Averaging over 200 MD trajectories 3 10"^ t.u. AT = O.OIT, 7 = 1. 
Dashed line - harmonic approximation at T_ — )• 0. 

linear in the limit AT — )■ when — )■ 0. Thus any quadratic mean values should be of the 
order ~ AT. 

Two temperatures of the background process T_ were tested: Ti = 0.2 and T2 = 5 
(from here we omit subindex '-' for brevity). Mean value ([xj^(t)]^) was analyzed as an 
example and results are shown in Fig. H] (fully identical properties have all quadratic values 
(correlators) of the types ( (t)x](t)) , ( Xj^(t)x](t)) ( Xj^(t) x](t))). As one expects, the 
quadratic form ([xj;(t)]^) linearly depends on AT: ([a;i(t)]^) ~ AT at Ti = 0.2. But the 
case is quite different at T2 = 5: mean value ( [xj;(t)]^) tends to a stationary value 0.064 in 
the limit AT — i- 0. It means that there exists some undamped stationary process x^(t) at 
high temperatures T even in the limit AT — ?■ 0. These results also can imply an existence 
of a threshold temperature Tthr separating two regimes - damped at low temperatures and 
undamped at high temperatures. 
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FIG. 4: Dependence of the mean value of quadratic form ((a;|)^) on the temperature difference 
AT. Fihed circles: T = 0.2 , empty circles: T = 5. Asymptotic value {(^i)'^) ^x^o ~ ^■^^'^ 
T = 5 (coefficient of linear regression 0.9993). Averaged over 100 MD trajectories 10^ t.u. each. 
= 5. The range of AT: 10"" < AT < 2-10-1 

IV. THE THRESHOLD TEMPERATURE 

Any process ^^{t) damps out at low temperatures and flattens out to a stationary value 
at higher temperatures even in the limit AT — )■ 0, and the temperature T of process x°(t) 
determines the different damping rates. And an illustrative process xi(t) was analyzed to 
determine an existence of a threshold temperature and its value ('tilde' marks the process 
xi(t) at AT = to avoid confusions). 

Process x^(t) can be exited in some or other manner. Usually x} and vj get random 
values in such a way that — 0)]^ — \ XlJ'^iH^ = 0)]^ = The particular choice 

of initial conditions does not influence the flnal results. 

Stochastic differential equation for the process x^(t) are 



dU _ dU^ 

dxj dxi 



— ^iix^ — 6ii\fXj^; 
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(9) 
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FIG. 5: a) Exponential damping of process x^(t) at different temperatures: T = 3.5 (circles), 
T = 3.8 (squares), T = 4.0 (triangles up), T = 4.2 (triangles down). Solid lines - linear regression. 
Averaging time ~5 000 — 10 000 t.u. 20 trajectories x'' were used to estimate the standard error, 
b) Damping coefficient {—a) as the function of temperature T of process x'^(t). Damping stops 
(a = 0) at Tthr ^ 4.07. iV = 5. 

with random forces determined by the difference of processes x(t) and x°(t), and damping 
at the extreme left and right oscillators; U and f/° are potential energies with coordinates 
x(t) and x°(t), correspondingly. Stochastic dynamics is implicitly ruled out by the 
temperature T of process x°(t). 



A. Two methods to find Tthr 

We consider the case of small temperature T when process x^(t) is damped out. The 
damping is determined by the viscous friction of left {—Xi) and right (— x^) oscillators in 
iQ. Gradually increasing the temperature we find its threshold value when process x^(t) 
becomes undamped. 

The damping of mean squared displacement of the first oscillator x\(t) was calculated. 
It was found that this process exponentially decays ( [xi(i)]^) oc exp(— at) and a depends 
on T (see Fig. [5^). One can see that the damping stops in the range 4.0 < T < 4.2. The 
dependence of coefficient a on the temperature T of process x^ (t) is shown in Fig. [5]d and 
Tthr ^ 4.07 at a = 0. 
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FIG. 6: Stationary values ((xj)^) at temperatures higher then Tthr in log-hnear coordinates. 
Time of averaging 10^ t.u. The temperature dependence was approximated by the function 
([xi(t)]2) ~ exp[-6/(r - Tthr)] at r > Tthr (sohd hne). Tthr ^ 4.09. iV = 5. 

Now we find the threshold temperature going "from up to down", going from higher 
temperatures. At high temperatures there exists the stationary process outcoming from 
random forces $ (see ([7]) - expression in square brackets). Process x^(t) decreases in the 
sense that all quadratic mean values tend to zero as temperatures approaches Tthr- When 
the threshold temperature reaches it threshold value, process x^(t) disappears (see Fig. |6]). 
The found threshold temperature is Tthr — 4.09. 

B. Time-resolved dynamics of process x^(t) 

To elucidate the reasons of strange dynamics of process x^(t) at high temperatures T 
we analyzed it more thoroughly. As above, A(t) = [5;J(t)]^ was calculated but without 
averaging over time. Results are shown in Fig. [7] for three temperatures T of process x°. 
One can see that A(t) behaves highly irregular. And numbers and heights of observed peaks 
increases with the growth of temperature until it becomes chaotic at high T. The mean 
values (A(T))^ at different temperatures T averaged over time interval r = 10^ t.u. increase 
with temperature. Mean values (A(t)), shown in horizontal solid lines, are nothing else 
than the stationary values calculated above. It was specially checked out that the dynamics 
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FIG. 7: Dependence of A(t) versus time at different temperatures of process x'^(t). a) T = 3.9; 
b) T = 4.3; c) r = 7.0. = 5, integration step h = 0.01. Oscillations at small times (insert 'A' to 
panel a) decay on average exponentially. Detailed shape of "excitation" at t w 93 800 is shown in 
insert 'B'. Mean values {'^{t))f^'Z^ are shown in horizontal solid lines. 

observed in Fig. [7] is not due to numerical artifacts. 

C. Heat conductance at small temperature gradients 

Our main concern is the computation of heat conductance at small temperature gradients. 
With this in mind we analyze an expression for the heat current in more details. And this 
analysis can also shed some light upon the problem why process x^(t) behaves in such strange 
manner. Remind an expression for force in the heat current: = — — Xt) 

is the force acting on the {i + l)th oscillator from left to right, and the derivative of the 
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potential energy u{xi+i — Xi) between oscillators is taken with respect to Xj+i. Then the 
expression for the local heat current ([8]) can be rewritten as 



Ji^i+i = ([Fi^i+i(x) - Fi^,;+i(x°)] xl^^) = 

{ — Xi) — (xi+i — XiY + {x^+i — x^) + — x^Y] xj+i) 



(10) 



where x}_^_i - velocity of {i + l)th oscillator in process x^. Difference of forces (in square 
brackets in the second line) is the polynomial of the third degree in the square root of 
temperature difference a/ AT (process x-*- is of the order of ^/AT as discussed above), and 



taking into account that the velocity xj^i is also of the order of V AT, it is the polynomial 
of the forth degree in a/ AT. But the coefficient of heat conductance is determined by the 
relation J/ AT therefor terms of the third and forth orders can be neglected at small values 
of AT. Then f lTOj) is simplified to 



J,;- 



i^i+i — ^ {xi+i a;^) 1 + 3 {x^^i Xj) Xj+i^ . (11) 

It is significant that the total heat current in f lTOj) at large temperature T and small 
temperature gradient AT is the difference of finite terms from processes x(t) and x°(t), 
vanishing in the limit AT — t- 0. And it is the reason why direct MD simulation is highly 
inefficient in this case and gives very large fiuctuations. But, as will be shown below, there 
exists an efficient method to overcome this difficulty. 

The behavior of process x^(t) is explained by the fact that it is determined not only 
by random Langevin forces ~ V AT, but also (and more essentially) by time- dependent 
random forces $j = [dU /dxi — dU^/dxi] (see ([7])). The plateau for the correlator ([xj(t)]^) 
equal to 0.064 at AT — i- is determined exclusively by random forces $j from stationary 
process x° (an illustrative example of one variable is considered in Appendix [B]). Thus, the 
dynamical process ^^(t) becomes the stationary one, determined by the background process 
x'^(t) at high temperatures. 



V. THRESHOLD TEMPERATURE IN THE LIMIT AT ^ 

In this section process x^ (t) is considered at an arbitrary temperature T and in the limit 
AT — !■ 0. Remind that process x^(t) ~ ~ V AT is completely suppressed at T < Tthr- To 
realize the limiting transition AT —t- in ([7]) it is convenient to divide both sides by V AT. 
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Then new coordinates are y{t) = x^{t)/\/ AT. It is also convenient to introduce normalized 
to unity random force 6 = / V AT. Then the linear equation for y{t) can be obtained 
as quadratic and cubic forces can be neglected (see ffTOj) ). The corresponding equation for 
y{t) can be derived if to rearrange one term in potential energy ([1]) keeping in mind that 
X, = xO + xi. ThenM(x,-Xi_i) = \[{x^i - xti) + {x\ - xU)f + \[{x^, - xti) + {x] - xU)]\ 
Transforming variables to y and retaining terms quadratic in y, one can get the potential 
energy in the form 

« = ^ E^^(^) - 3^{t) = 1 + 3 [x°(t) - xU{t)f , (12) 

i 

where gi{t) are time- dependent random coefficients of rigidity determined by the dynamical 
process x°(t). It is illuminating to note that the problem of heat conductance can be 
reduced to the quadratic potential energy in the limit AT 0. Corresponding SDEs have 
Langevin source with unit temperature at the left oscillator and zero temperature at the 
right oscillator: 

Vi = -9i{yi - Vi-i) + 9i+ii:yi+i - Vi) + ^ni^ - m) - ^iNVN ■ (13) 

It should be also noted that if the ID lattice with an arbitrary interaction (Morse, Toda, 
LJ, etc) is analyzed then the corresponding equation will be the same, and random rigidities 
are Qi = U"{x^ — x^^^) where U is some or other type of potential energy. Equation for 
arbitrary systems (with arbitrary neighbor radius of interaction) can be also written in the 
general form as 

M 

yi = -~Yl + i^-yi)- ^iNVN , (14) 

where A^^- - matrix of second derivatives of potential energy depending on x*^, and M is the 
number of neighbors. Equation (IT^ is valid for arbitrary systems. 

Equations f|T4|) define the stationary random process only if temperature T < T^hr- As 
temperature approaches the value Tthr, quadratic mean values diverge. It is shown in Fig. [HI 
And this is the third method to find Tthr- 

There was analyzed the case of low temperatures T when process x° is "weak". Then 
rigidity coefficients gi are close to unity (see (fT2|) ). And as an example we consider the 
lattice where actual rigidity coefficients Qi f|T2l) are substituted by the mean value taken 
from the equilibrium Gibbs distribution = Qq^T) and go{T) = 1 + 3 ((x° — x°_i)^). This 
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FIG. 8: Dependence of mean value ((yi)^) vs. temperature T. At temperature T close to Tthr ~ 4.1 
process diverges. Averaging over 20 MD trajectories 210^ t.u. length each, ((yi)'^) ~ 0.61 at 

harmonic model is exactly solvable and results are shown in Fig. |9]0ne can see that process 
y{t) is damped out in the model with constant rigidity in contrast to the case when actual 
values f|T2|) are used. And one can conclude that the growth of process y{t), when temper- 
ature increases, is determined by an increase of fluctuations but not only by the increase of 
rigidities. 

Process y{t) diverges at high temperatures. And it gives one more possibility, the fourth 
one, to find the threshold temperature. To attain this end the equilibrium process x''(t) 
at temperature T is established. Then process y(t = 0) is excited in some or other way 
(its initial conditions do not influence the final results). And the evolution of the y{t) is 
analyzed. One can see (Fig. [10]) that the process exponentially damps out at T < Tthr and 
exponentially grows at T > Tthr- 

It should be stressed out that the method just described differs from the previous one 
(Fig. [6] and discussion) . The nonlinear case was considered there and its stationarity was 
conditioned by nonlinear terms in forces which are absent in the harmonic approximation. 

Four methods give the threshold temperature Tthr ~ 4.1. This temperature was found 
for a fixed lattice length = 5. Larger lattice lengths were considered and the dependence 
of Tthr on the lattice length A^ is shown in Fig. [TTl Approximate fitting gives dependence 
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FIG. 9: Dependence of the mean squared displacement ([yi(i)]^) vs. temperature T. Circles - 
MD simulation of SDEs ()13p ; solid line - model of mean rigidities in the harmonic approximation. 
Averaging over 20 trajectories 2 10^ t.u. each. 




FIG. 10: Exponential dependence of the quadratic form oc exp(— at) on time. Tempera- 

tures from bottom to top: T = 3.5, 3.7, 3.9, 4.1, 4.3, 4.5. Averaging over 20 MD trajectories, 2 10^ 
t.u. each. The dependence of coefficient a on temperature is shown in insert. T^^r ~ 4.1 is found 
from the condition when a = 0. 



17 



-3 -I 1 1 1 1 1 1 

0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 

ig N 

FIG. 11: Dependence of Tthr vs. lattice length N in log-log coordinates. Solid line is the linear 
fitting Tthr ~ N~^. 

Tthr ~ 6 ■ 10^ A^"'^. It means that the majority of usually studied lattices are in the state 
when their temperatures are much higher then the threshold temperature (e.g. if > 100 
then Tthr < 6 ■ 10"^). 

VI. SOUND VELOCITY AND SOLITONS IN /3-FPU LATTICE 

Heat conductance is observed in both regimes, - higher and below the threshold temper- 
ature. And an attempt was undertaken to find the soliton contribution to the heat conduc- 
tance at T > Tthr- With this in mind, the correlator {Axk{t) Axk+m{t + ''")) was analyzed 
(Axj(t) is the displacement of ith oscillator from equilibrium at time instant t). In numerical 
simulations we fixed the time shift r = 20 t.u. and calculated the corresponding correlator 
(A^ = 101, T = 2). Results are shown in Fig. [121 The correlator {Ax^oit) Ax^o+mit + 20)) 
has peaks at the coordinate shifts m = ±25. It allows to calculate the velocity of excitation 
propagation t>exc and v^xc ~ 1.25. This velocity is higher then the sound velocity calculated 
in the harmonic approximation v^^^^^ = 1 at /3 = 1. 

Initially these peaks were attributed to solitons. But more thorough analysis shows that 
this concepts is not valid. Let we have the /3-FPU potential u{y) = |?/^ + jy"^ {(3 = 1 and 
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FIG. 12: Correlator (Ax5o(t) Ax5o+m(i + 20)) versus the lattice position (number). /3-FPU lattice 
of = 101 oscillators is used. T = 2. Arrow shows the site position n = 50 



Hi = Xi — Xj_i). In [SQ!] it was shown that there exists a spectrum of frequencies which are 
proportional to the harmonic ones, according to a well defined law. Therefor the /3-FPU 
potential can be represented as 



and an expression in brackets can be replaced by an effective harmonic rigidity 

u{y) = hs^y'^ 



(15) 



(16) 



The problem is to find k^s- It can be done in terms of a mean field approximation (MFA). 
Mean value of potential energy is 



(17) 



where (y^) is the mean value of y^. 

In the harmonic approximation (at not too high temperatures) mean values of potential 
and kinetic energies are equal (up) = (wk)- In canonical ensemble the identity (wk) = T/2 is 
valid for ID systems. Then 

{y') = I (18) 
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FIG. 13: The dependence of VqS versus temperature T: solid line - dependence ([2^ : empty circles 
- MD simulation. 

The self consistency of the MFA is (expression in brackets in f|T5|) = keg) 

1 + 1 {y')) = (19) 

From (fT8|) it follows that (y^) = T/k^s and substitution of (y^) = T/k^.^ into (fT9!) gives the 
self-consistent equation for kes 

1 + r/(2fcefr) = k,s (20) 

with the solution 

1 /l T 

Thereby the the effective ( "nonlinear" ) sound velocity 




Ves = Vhs = \/ o + \/ 7 + 77' ("^ = 1) (22) 



is higher then the velocity f ^ouJJJj = 1 found in the harmonic approximation. Eq. [22] gives 
Ves = 1-27 for T = 2 what coincides with the value fexc ~ 1-25 found from correlation 
functions. The dependence of effective sound velocity versus temperature is shown in Fig. [T3i 
Note that the MFA is valid up to very high temperatures T = 10, while this approach 
originally is well suited only for low temperatures, and the effective sound velocity exceeds 
its harmonic value (at T = 0) by > 50%. 
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FIG. 14: Solitons and breasers running out of the lattice. A - breather, B - sohton, C - pair 
of antisohtons (solitons of elongation). Arrow at n = 200 shows the border separating initially 
thermalized and "cold" parts of lattice. Initial temperature of the left part (1 < n < 200) of the 
lattice T = 10. 

Next we try to find direct evidences on the solitons participation in energy transfer. It 
was done in the following manner. Initially lattice of = 200 oscillators was thermalized for 
some time to reach the thermodynamic equilibrium. Then the "cold" lattice (with zero ve- 
locities and displacements) with 1000 oscillators was switched to the right end of the lattice. 
Solitons, if they exist in the initial lattice, should "run out" to the cold lattice. The same 
is valid for the moving breathers. (Note that in the continuum approximation the mKdV 
equation corresponds to the discrete /3-FPU potential. And one can find analytical expres- 
sions for solitons of compression, antisohtons of elogation and different types of breathers 



m 



60| ). We waited some time till excitations run out of the lattice to its cold part where 



they can be observed. Results are shown in Fig. [H] 
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Analogous approach to visualize breathers in 2D lattice with three different on-site po- 



tentials was utilized in 



61| where initially thermalized lattice was cooled from the borders 



and breathers were detected after thermal noise was deleted through damping boundaries. 



The possibility of energy transfer due to solitons was conjectured three decades ago 



3Q|. 



Less studied is the possibility of energy transfer by breathers. One suggested mechanism 



is the Targeted Energy Transfer 6^, |63|] when an efficient energy transfer can occur under 



a precise condition of nonlinear resonance between discrete breathers. Various aspects and 
possible applications of energy transfer by breathers are considered in 64l |. 



VII. CONCLUSIONS 



In conclusion we briefly summarize our main results. A new method is developed which 
allows considerable decreasing of the computation time in calculations of the heat conduc- 
tance at low temperature gradients (temperature T + AT of the left lattice end and T - of 
the right end and AT/T ^ 1). This success was achieved by the separation of the total 
process of heat conductance into two parts: an equilibrium process x°(t) at equal temper- 
atures T of both lattice ends and non-equilibrium process x^(t), responsible for the energy 
transport, which occurs at temperature AT of one end and zero temperature of other end. 
The equilibrium (background) process strongly influences the transport properties: there 
exists the threshold temperature Tthr above which some undamped characteristics are ob- 
served; more precisely, correlators of the types (^x}{t)x^{t)') do not tend to zero at AT — )■ 0, 
as expected, but have certain nonzero values; at T < Tthr "normal" dependence is observed, 
i.e. these correlators have zero values when AT 0. The reason of two distinct behaviors is 
not due to the temperature of the background process x°(t) but sooner to the temperature 
fluctuations. An illustrative example of one variable is briefly analyzed where the threshold 
temperature is also found. The model of one variable has a rich family of solutions depending 
on the parameters and serves to be investigated in more thoroughly. 

The threshold temperature was found by few methods and scales ~ A^~^ with the lattice 
size A^. All practically interesting systems lies above Tthr- The threshold temperature is 
not sharply pronounced and arbitrarily separates two mechanisms of the heat conduction: 
the phonon mechanism prevails at T < Tthr, and at T > Tthr the soliton contribution starts 
to play more significant role with the increase of temperature. Highly probable that the 



22 



temperature fluctuations are responsible for the solitons generation. 

Analytical solutions for solitons and breathers are known for the ^-FPU lattice and these 
excitations were directly observed in numerical experiments. Our findings on the soliton 
contribution to the heat conductance are in accordance with the general scenario of heat 
conductance: phonons gave main contribution to the heat conductance at low temperatures 
and solitons more and more dominate when temperature increases. 

We found no relations between the well known weak and strong stochasticity thresholds 
and the threshold temperature in the present paper: they have different energy ranges and 
different dependencies on A^. Additional difference is due to different statistical ensembles 
used: traditionally stochasticity thresholds are found in microcanonical ensemble, but criti- 
cal temperature is observed in canonical ensemble where energy equipartition is realized at 
any temperature. Statistical properties do coincide in the thermodynamical limit N ^ oo 
for //-canonical and canonical ensembles, but the dynamical properties can differ. 
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FIG. 15: Standard errors of the heat current computed by solving SDEs ([2]) (empty circles) and 
new SDEs (filled circles). Temperature T = 0.1, AT = O.OIT. Averaging over M = 100 

MD trajectories 10^ t.u. each; total time of computation 10^ t.u. 

Appendix A: Comparison of efficiencies of diff"erent methods in computation of heat 
conductance 

The main problem in the calculation of heat conductance is to find the heat current with 
an appropriate standard error. And we consider below the effectiveness of the separation 
of the total process into the sum of two: x(t) = x°(t) + x^(t) in the sense of computer 
time expenditure (see ([2]) and (EI)- ([7]) ), where x° is the equihbrium background process and 
process x^ is responsible for the heat transport. 

The comparative efficiency of two approaches ('old' - solving of SDEs ([2]) and 'new' - 
two SDEs ([6])- ([7]) ) to the calculation of the heat flux can be estimated as the relation of 
their standard errors 6 at equal conditions of computation: Eff = ^oid/f^ncw The result is 
shown in Fig. [151 One can see that the standard error is systematically less in the suggested 
approach as compared to the usually utilized. 

More impressive is the behavior of efficiency at decreasing of the temperature gradient AT 
(see Fig. [T6|) and at equal parameters of numerical simulations. Here it should be emphasized 
that standard errors increase with the diminishing of AT in wide range 10~^ < AT/T < 10~^ 
and the efficiency Eff > 1 remains as before. The standard error 6 increases as 5 ~ T/AT. 
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FIG. 16: Standard errors vs. inverse temperature gradient AT. Filled circles ~ results of solution 
SDEs ^ and filled circles - solution of ©-([T]). T = 1, Af = 100 MD trajectories 10"^ t.u. length 
each in both cases. 

Appendix B: An example of one variable 

We consider an example of an equation with one variable x for the harmonic oscillator 
with damping 

X = —k(t)x — ■jx . (Bl) 

where k{t) is the stochastic rigidity. This equation is the illustrative analogue of multi- 
variable SDEs equations ([7]) for the process ^^(t). The variable substitution xexp(— 7^/2) — )• 
X eliminates the damping and flBl|) can be reduced to 

X = -k{t)X. (B2) 

Potential energy f|T2l) has the form u = \g{t)y'^ in the case of one variable, where g{t) = 
1 + 3 x^(^) S'lid x{t) - stochastic process generated by the background process x°(t). And it 
is reasonable to choose the random rigidity in (IB2p in the form 

k{t) = l+ ez^{t) , (B3) 

where e is free parameter and z{t) is the stationary random process describing the dynamic 
of the harmonic oscillator influenced by Langevin source with temperature T: 

z = -z + ^ - -fz (B4) 
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and((e(ti)e(t2) = 27T5(ti-t2)) 

We compare the solution X{t) flBll) - flB3P with the solution of well known deterministic 
Mathieu equation 

y = -[l + gcos\t)]y. (B5) 

Different types of solutions of the Mathieu equation depend on the parameter g and initial 
conditions. There exists such Qcr that the solution is the sum of periodic functions a.t g < gcv, 
and the solution is the superposition of periodic functions multiplied by the exponentially 
increasing and decreasing function exp(±/it) at g > gcr- 

Equations (lB2p - (lB4p also have a rich family of solutions depending on initial conditions 
and parameter values. As an illustrative example we consider the following set of parameters: 
X{t = 0) = 0.5, X{t = 0) = 0, e = 50, 7 = 1. 

Below we demonstrate only the qualitative behavior of process X{t) depending on the 
temperature T of stochastic process flB4|) and results are shown in Figs. [TTT a-c). One can 
see different regimes as T increases. And there exists some critical temperature Tcr above 
which the process X{t) diverges (Tcr ~ 100). It should be noted that the overall scenario 
strongly depends on the choice of initial conditions {X{t = 0), X{t = 0)) and the particular 
sequence of random Langevin forces {^} in (1B4I) . At larger times process X{t) becomes more 
complex. The full analysis of system flB2|) - flB4|) is not our primary goal, but these equations 
serve more intensive attention. 
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